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Abstract 

Reported here are resuits of theoreticai calculations on Fe(H 2 0)6 3+ , 
Fe(H 2 0) 5 (OH) 2 +, three isomers of Fe(H 2 0) 4 (OH) 2 + , and Fe(H 2 0) 3 (OH) 2 + 
to investigate the molecular mechanisms of hydrolysis of ferric ion in wa- 
ter. The combination of density functional electronic structure techniques 
and a dielectric continuum model for electrostatic solvation applied to the 
Fe(H 2 0)6 3+ complex yields an estimate of -1020 kcal/mol (experimental val- 
ues -1037 kcal/mol to -1019 kcal/mol) for the absolute free energy of the 
aqueous ferric ion. The free energy change for the first hydrolysis reaction is 
predicted reasonably (2 kcal/mol predicted compared to 3 kcal/mol experi- 
mental). For the second hydrolysis reaction, we found an unexpected low en- 
ergy isomer of Fe(H 2 0)4(OH) 2 + with five ligands in the inner sphere and one 
water outside. The hexa-coordinate cis and trans isomers are, respectively, 
slightly lower and higher in energy. Calculations on the penta-coordinate 
species Fe(H 2 0)3(OH) 2 + suggest that extrusion of the outer sphere water is 
nearly thermoneutral. The reaction free energy for the second hydrolysis is 
predicted in the range 16-18 kcal/mol, higher than the experimental value of 
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5 kcal/mol. In view of the facts that the theoretical predictions are higher 
than experimental values, and that novel structures were encountered among 
products of the second hydrolysis, we argue that conformational entropy is an 
important omission in this theoretical treatment of net reaction free energies. 
A fuller cataloging of low energy hydrolysis products and direct calculations 
of partition functions of the isolated complexes should help in modeling equi- 
librium speciation in groundwaters. 
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1. Introduction 

A molecular theoretical account of the free energetics of the reactions 

[Fe(H 2 0) n (OH) (m _ 1) ] 4 - m + H 2 -> [Fe(H 2 0) n _!(OH) m ] 3 - m + H 3 0+ (1) 

in water can benefit from discrimination of the conformers that are present under common 
conditions. Entropic contributions to the solution thermodynamics may reflect multiple 
configurations that occur. Thus information on the conformers present should assist in 
accurately describing temperature variations of solvation properties. In addition, theoretical 
study of the molecular structure of the species participating these reactions should teach us 
about the molecular mechanisms involved and provide a benchmark of current theoretical 
tools for modeling speciation of metal ions in groundwaters • 

The condition of ferrous and ferric ions in solution has long been of specific interest to 
explanations of electron exchange processes @J^,|,0. The hydrolysis product Fe(H 2 0) 5 OH 2+ 
often contributes significantly to the rate of ferric ion reduction in water because the specific 
rate constant for ferric-ferrous electron exchange is about a thousand times larger for the 
hydrolyzed compared to unhydrolyzed hexaaquoferric ion ||. The net standard free energy 
change (about 3 kcal/mol |J) for the first hydrolysis reaction, 

Fe(H 2 0) 6 3+ + H 2 - Fe(H 2 0) 5 (OH) 2 + + H 3 0+, (2) 

is small compared to the size of the various contributions that must be considered in theo- 
retical modeling of these solution species. 

In contrast to the great volume of work on electron transfer involving such species, 
reactions of particular interest here transfer a proton from a water ligated to a Fe 3+ ion to 
a free water molecule. The second hydrolysis, 

Fe(H 2 0) 5 OH 2 + + H 2 -> Fe(H 2 0) 4 (OH) 2 + + H 3 0+, (3) 

raises more seriously the possibility of participation by several isomers in the mechanism 
and thermodynamics of these reactions. 
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This work takes ferric ion hydrolysis as a demonstration system for determination of 
the utility of current theoretical tools for modeling of speciation of metal ions in aqueous 
solution. At this stage we have not considered multiple metal center aggregates. We report 
below encouraging agreement with experimental thermochemistries but with unexpected 
results that must be addressed for further progress in describing these chemical systems in 
molecular terms. 

2. Approach 

The physical idea of treating an inner solvation shell on a different footing from the rest 
of the solvent has considerable precedent [ID . Friedman and Krishnan referred to these ap- 



proaches as hybrid models and criticized them on the grounds that individual contributions 
to the thermodynamic properties of final interest could not be determined with sufficient 
accuracy to draw new conclusions from the final results. However, particularly for highly 
charged and chemically complex ions, the hybrid approaches seem indispensable |TT|. Com- 
putational and conceptual progress in the 25 years since that Friedman-Krishnan review has 
made these approaches interesting once again. In what follows, we first elaborate on the 
calculational techniques used. An appendix discusses some of the statistical thermodynamic 
issues underlying these hybrid models. 

2.1. Electronic Structure Calculations 

Reliable estimation of free energies of dissociation in aqueous media requires a similarly 
reliable determination of the process in the gas phase. Given the well known difficulties 
associated with electron correlation in transition metal complexes, and our intent to extend 
this work to larger and more complex systems, we chose to explore this problem with the 
B3LYP hybrid density functional theory (DFT) [fT^JTcjl . This approximation is an effective 
compromise between accuracy and computational expense, and has been shown to give 



usefully accurate predictions of metal-ligand bond energies in a number of molecules |1J,|15 



In the present application however, we were less interested in the metal-ligand bond energy 



than in the energies for successive elimination of protons from the H2O ligands bound to 
Fe 3+ . This may be a somewhat less severe test of the functional. The details of our approach 
and estimates of the errors follow. 

The geometry and vibrational frequencies of the metal complex were determined using 
the 6-31 1+G basis set for the metal ion and the 6-31G* basis set for the ligands [|13||. The 



former is contracted from Wachters' |T6[ primitive Gaussian basis (including the functions 



describing the 4s and 4p atomic orbitals) and augmented with the diffuse d-function of Hay 



17| . The ligand basis set contains a polarization function of d-character on the oxygen. At 
the optimum geometry, a single calculation was done to determine the energy in a more 
extended basis (6-31++G**) that includes polarization and diffuse functions on both the 
oxygen and hydrogen centers. Atomic charges were determined in this extended basis as 
well using the ChelpG |18| capability (Rp e =2.02 A ) in the Gaussian94 package [13]. All Fe 



species are high-spin (d 5 ) treated in the spin-unrestricted formalism. 

One might argue, on the basis of the OH - character expected in the product complexes, 
that the geometry optimization should also be done with a diffuse basis set for the ligands. 
This was found unnecessary by explicit calculations on the partial reaction 

Fe(H 2 0) 6 3 + -> Fe(H 2 0) 5 (OH) 2 + + H+. (4) 

The endothermicity of this reaction computed from a single point calculation with the 6- 
31+G* basis at its optimum geometry differs by only 0.2 kcal/mole from the result obtained 
with the smaller basis. More important are additional polarization and diffuse functions on 
the hydrogens; for example, neglecting zero-point corrections, the B3LYP endothermicities 
are 29.2 kcal/mole(6-31G*), 25.4 kcal/mole(6-31+G*), and 28.1 kcal/mole(6-31++G**). 

A final consideration regarding basis set convergence is the importance of polariza- 
tion (f-functions) on the metal center. This was examined by augmenting the 6-311+G 
Fe basis with two f-functions (a = 0.25,0.75) split about a value (a = 0.4) optimized for 
Fe(H 2 0)6 3+ . In this larger basis an endothermicity of 27.9 kcal/mole results compared to 
the value 28.1 kcal/mole obtained without the f-functions. The correction to the second 



hydrolysis is only slightly larger, ~ 0.5 kcal/mole. We thus conclude that the B3LYP de- 
protonation energies are nearly converged with respect to further improvements in the basis 
set. A conservative estimate might be associated with the change observed between the 6- 
311+G/6-31+G* results and the 6-311+G(2f)/6-31++G** results, or about 2-3 kcal/mole. 

Convergence of the results with respect to basis set does not address the accuracy ex- 
pected with the B3LYP method. In so far as the water molecules retain their identity when 
bound to the metal ion, we expect the error in the hydrolysis reaction to be approximated 
by the error observed for the analogous reaction for a single H2O molecule. The procedure 
outlined above was therefore used to determine the energies and zero-point energies (ZPE) 
for a number of species associated with the neutral and ionic dissociation channels of H 2 0. 
They are reported in Table I, and can be used to determine the properties and reaction 
energies reported in Table II. The H 2 thermochemistries shown there are in excellent 
agreement with experiment. For example, the error associated with the reaction 

H 2 -»• H+ + OH" (5) 

is underestimated by less than 2 kcal/mole. The errors in the other dissociation channels 
are similar. 

2.2. Electrostatic Solvation Calculations 

Electrostatic interactions of hydrated ferric ions with the aqueous environment are ex- 
pected to be of first importance in this solution chemistry. These hydrolysis reactions are 
written so that other contributions to the net free energy, e.g. packing, might balance 
reasonably. Thus, continuum dielectric solvation modeling is a reasonable approach to 
studying these solution processes; it is physical, computationally feasible, and can pro- 
vide a basis for more molecular theory ■ The dielectric 



model produces an approximation to the interaction part of the chemical potential, fi* of 
the solute; in notation used below /1* = — RT ln(exp(— AU))o where for this case AU is 
the solute-solvent electrostatic interaction potential energy and the brackets indicate the 



average over the thermal motion of the solvent uninfluenced by the charge distribution of 



the solute |T9| , pO| , pT1 , p2| , p3| , pl| , p7| , p9| . Free energy contributions due to non-electrostatic in- 
teractions are neglected. We address free energy changes due to atomic motions internal to 
the complexes on the basis of gas phase determinations of the vibrational frequencies and 
partition functions. The solvation calculation treats all species as rigid; we comment later on 
the consequences of this approximation and on the possibilities for relaxing it . The solute 
molecular surface was the boundary between volume enclosed by spheres centered on all 
atoms and the exterior. The sphere radii were those determined empirically by Stefanovich 
pOfl , except Rp e =2.08 A for the ferric ion. The ferric ion is well buried by the ligands of 
these complexes and slight variations of this latter value were found to be unimportant. 
The numerical solutions of the Poisson equation were produced with a boundary integral 



procedure as sketched in references pqj27| typically employing approximately 16K boundary 
elements. The accuracy of the numerical solution of the dielectric model is expected to be 
better than a kcal/mol for the electrostatic solvation free energy. 

3. Results and Discussion 

The pertinent energies for Fe(H20)6 3+ and other species involved in hydrolysis processes 
are compiled in Tables I and II. The magnitudes of the various components entering into the 
gas phase hydrolysis free energy are reported in Table III, while Table IV adds the solvation 
contributions to the free energy and compares the total with experiment. Table V summa- 
rizes geometrical information and partial atomic charges computed from the electrostatic 
potential fit are presented in Table [VT]. 



3.1. Fe(H 2 Q) 



3+ 



The structure of the Fe(H20)6 3+ complex, T^ symmetry, is shown in Figure |T[ The 
B3LYP approximation gives an Fe-0 distance of 2.061 A . By way of comparison, with 
the identical basis set the Hartree-Fock approximation yields 2.066 A . These distances 
are similar to the recent Hartree-Fock results of Akesson, et al. ||31||,Rjr e o = 2.062 A , 



and the gradient-corrected DFT(BPW86) calculations of Li, et al. |32| , R^ e o= 2.067 A . 
These theory values cluster in the upper end of the range of distances (1.97 A - 2.05 A ) 
determined experimentally p3, 34|,35| . Neutron scattering measurements report 2.01 A in 
concentrated electrolyte solutions [33]. The solution EXAFS result (1.98 A ) is close to the 



crystallographic determinations, 1.97 A - 2.00 A [[33 ]. In their recent review, Ohtaki and 
Radnai discuss a number of other experiments and conclude the distance lies in the range 
of 2.01 - 2.05 A []35|. 

Though these theory values are compatible with the upper end of this range, they are far 
from the EXAFS result (1.98 A ) that might be the most reliable experimental determination. 
While investigating this question, we found that the B3LYP value is stable with respect to 
further improvements in the basis; if the basis is augmented with an f-function on the metal 
and diffuse functions on the oxygens, the equilibrium distance decreases only slightly {Kfs-o 
= 2.053 A ). As noted above, the B3LYP result is in close agreement with that from the 
Hartree-Fock approximation. It would not be surprising, however, if the HF approximation 
overestimated the bond length. Akesson, et al. p6| , p7| have previously observed that the 
HF bond lengths for a series of first row transition metal hydrates are systematically too 
long. They find that correlation effects and the influence of the second hydration shell each 
act to reduce the bond length by about 0.01 A . Even with these corrections the theoretical 
bond length is still much larger than 1.98 A . As an additional data point, the local- density- 
approximation with the present basis yields RFe-o — 2.013 A , in better agreement with 
the EXAFS and crystallographic determinations. However, while the LDA generally gives 
reasonable geometries for transition metal complexes []5§| , |3"5[] , Sosa, et al. J55| find that it 
tends to underestimate the lengths of dative bonds, such as the one discussed here, and this 
might further support a bond length toward the upper end of the range. 

The atomic partial charges computed here for the hexaaquoferric ion are given in Table 



VI. Note that the charge assigned to the ferric ion is substantially less than 3e, and that 



the magnitudes of the charges on the oxygen and hydrogen atoms are substantially greater 
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than is common with effective force fields used for simulation of liquid water. 

As discussed in the appendix, a value for the absolute free energy of the ferric ion can 
be obtained from a free energy of formation of the Fe(H20)6 3+ complex from an isolated 
Fe 3+ ( 6 S) ion and six water molecules (Table III), provided that no solvation contribution is 
included for the atomic ion and that the actual, not standard state, species concentrations are 
used. That absolute free energy in aqueous solution is -1020 kcal/mol. Experimental values 
range from -1037 kcal/mol to -1019 kcal/mol |40|j4l|Jl0|j42| . We note that this theoretical 



value takes no account of solvation effects due to non-electrostatic interactions. 

This agreement with experiment is encouraging. In computing the absolute free energy 
there are a number of large contributions to the final result. This contrasts to the hydrolysis 
reactions where we have arranged things to encourage cancellation of errors. It is interesting 
to examine the components contributing to the free energy in the gas phase (Table III). Not 
surprisingly, the zero-point correction, +15 kcal/mole, is significant. Under the standard 
conditions, hypothetical p=l atm ideal gas with T=298.15 K, there is the large, unfavorable 
differential entropy contribution, +47 kcal/mole. This arises from the necessity of sequester- 
ing six water molecules in a dilute gas. The net free energy found for the gas phase reaction 
is -602 kcal/mole. Table IV addresses the solution phase aspects of this thermochemistry. 
(See the Appendix also.) About half of the ideal entropy penalty mentioned is regained in 
the liquid because of the higher concentration of water molecules. The differential solvation 
free energy adds another -391 kcal/mole of (favorable) free energy for a net absolute free 
energy of hydration of -1020 kcal/mole. 

This estimate of the absolute free energy of the ferric ion is close to the value - 



1037 kcal/mole reported by Li, et al. |32] as a hydration enthalpy. The Li, et al, value con- 
tains some terms that would be appropriate if the hydration enthalpy were sought and some 
terms that would contribute to the hydration free energy; thus comparison of that previous 
value with the present result is not straightforward. In fact, the most pragmatic calculation 
of the hydration enthalpy within the dielectric continuum model is nontrivial. The enthalpy 
would be then obtained by determining a temperature derivative of the chemical potential. 



The appropriate temperature derivative determines the enthalpy directly or, alternatively, 
the solvation entropy so that the desired enthalpy might be determined by differencing with 
the already known chemical potential. That temperature derivative would generally involve 



also the temperature variation of the radii-parameters |T9| , pO| , pT| , p7| , p9| . The variations with 
thermodynamic state of the radii-parameters have not been well studied |21[| . There is good 
agreement, however, in many of the components contributing to this hydration enthalpy. 
Li, et al., report a gas phase energy of formation of -652.2 kcal/mole with the BPW86 
gradient-corrected functional compared to the present B3LYP result of -655.2 kcal/mole, 
and their solvation free energy is -444 kcal/mole vs. the -441 kcal/mole found here. Such 
close agreement is encouraging, although somewhat fortuitous. Our gas phase energy of 
formation includes a correction, not considered by Li, et al., of some 15 kcal/mole for the 
zero-point energies. It is encouraging, however, that the estimates disagree by only about 
15 kcal/mole, or some 1.5%. Additionally, the definition of the molecular surface adopted 
by Li, et al., for the solvation calculation was substantially different from that used here, as 
were the radii-parameters used. 

3.2. First and second hydrolysis reactions 

We turn now to the first deprotonation reaction Eq. |2[ The structure found for 
Fe(H 2 0)50H 2+ is displayed in Figure 0. The Fe-0 (hydroxide) distance is 1.76 A and 
the Fe-0 (water) distances lengthen to 2.10 A - 2.15 A . Assembling the results for the 
standard free energy of this reaction we find 2 kcal/mol, in surprisingly good agreement with 
the experimental value of 3 kcal/mol ||. This computed net free energy change is composed 
of approximately -148 kcal/mol exothermic (favorable) change in isolated molecule free en- 
ergy and +150 kcal/mol (unfavorable) net increase in solvation free energy. The solvation 
contribution favors the reactant side here because it presents the most highly charged ion. 
Changes in the radius assigned to the iron atom in the range 2.06 A <Rp e <2.10 A lead 
to changes of ± 1 kcal/mol in the predicted reaction free energy. This reaction was also 
considered by Li, et al., J32] who find it to be exothermic by 14 kcal/mole. 
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To treat the next hydrolysis Eq. |3] we must consider Fe(H 2 0)4(OH) 2 + species. Figure 
3 shows the stable structures found. Further lengthening of both the Fe-0 (water) and 
Fe-0 (hydroxide) distances is noted(Table V) in the cis and trans six-coordinate species 
by 0.07-0.14 A compared to the Fe(H 2 0)50H 2+ . In the gas phase, the cis structure is 
predicted to be the lowest energy conformer, slightly (1 kcal/mole) below the trans isomer. 
This preference is reversed in solution, where the trans isomer is predicted to be slightly 
more stable. We note that the most current force fields applied to simulation of ferric ions 
in solution place the cis structure significantly higher in energy than the trans ||. 

We were surprised to discover a stable outer sphere complex during our search for the 
stable trans structure. The structure is given in Figure 3. The distances between the 
hydroxyl oxygens and the outer sphere water are typical of hydrogen bonds. Explicit calcu- 
lation of the vibrational frequencies show it to be a true local minimum, lying less than a 
kcal/mol higher in energy than the cis conformer. 

The interaction of the outer sphere water with the remainder of the ferric hydrate com- 
plex can be partially characterized by finding the energy of that complex without the outer 
sphere partner. The structure obtained for that penta-coordinate ferric ion is shown in Fig- 
ure 4. The penta-coordinate complex can stably adopt a conformation similar to that of 
Figure 3 also in the absence of the outer sphere water. In terms of the zero-point corrected 
electronic energy, this outer sphere complex is stable with respect to loss of the H 2 by 
about 7 kcal/mole. Consideration of the entropic contributions to the free energy find the 
complex still stable with respect to loss of H 2 0, but the differential solvation contributions 
reverse this conclusion. Thus, the dissociation of the outer sphere complex is an essen- 
tially thermoneutral process. We suspect such intermediates play an important role in the 
mechanism of the ligand exchange with the solvent. 

These three conformers lead to estimates of 16 kcal/mol, 16 kcal/mol, and 18 kcal/mol 
for the reaction free energy of the second hydrolysis for trans, outer, and cis sphere products, 
respectively; the experimental value is approximately 5 kcal/mol P,|4"3~|. We note that as in 



the case of the first hydrolysis, the gas phase predictions lead to an exothermic reaction; it 
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is the differential solvation free energy that tips the scales in favor of endothermicity. 



3.3. Role of conformational entropy 

The issue of the internal motions of the complexes, and the near degeneracy of several 
conformers, raises also the issue of conformational entropy of these species. For example, in 
the second hydrolysis reaction, factors such as -RTln3 arise if the three isomers discussed 
here are considered isoenergetic. [However, the multiplicity of isoenergetic states will surely 
not be just 3 and, furthermore, entropy differences are required here.] As another example, 
the Th hexaaquo complex surely has a number of low-lying structures involving the rotation 
of the plane of an individual H 2 0. More generally, this conformational entropy would be 
appropriately included by computing the solvation contribution to the chemical potential, 
fi*, of the complex according to 



fi* = -RT In xf>e-» 9 W ar (6) 

c 



where the sum indicated is over conformations c weighted by the normalized population, the 
mole fractions of conformers when there is no interaction between solute and solvent, 
and further the summand is the Boltzmann factor of the solvation free energies, perhaps 
from a physical model such as the dielectric model used here, for each conformation. [The 
treatment of the thermodynamics of flexible complexes in solution is also discussed further 
in the appendix and in Reference (44[].] The solvation contribution to the chemical poten- 



tial Eq. 10 would then be combined with the isolated cluster partition function to obtain 



the free energies of the species involved and free energy changes for the reactions ||44|| . Fi- 
nally, an entropic contribution, including any conformational entropy, would be obtained 
by temperature differentiation of the full chemical potential. The isolated cluster partition 
functions would properly include an entropy associated with the multiplicity of isoenergetic 
conformational states. It is reasonable to expect that this conformational entropy increases 
progressively with hydrolysis, i. e. the products here are "less ordered," and have higher con- 
formational entropy than their reactants. It is thus significant that the predicted reaction 
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free energies are higher than the experiment; inclusion of conformational entropy should 
lower these reaction free energies. 



We note that Li, et al. [g^] argue that inclusion of a second solvation shell substantially 
improves the accuracy of the thermochemical predictions. The present results do not seem to 
force us to larger clusters. Although it is true that proper inclusion of more water molecules 
should permit a more convincing treatment, inclusion of more distant water molecules makes 
the neglect of conformational entropy less tenable. In any case, comparison of dielectric 
model treatments with thermochemical results only permits limited conclusions because of 
the empirical adjustability of radii-parameters. 

4. Conclusions 

Given the balance of the large contributions that must be considered, the observed ac- 
curacy of the computed hydrolysis reaction free energies reactions is encouraging. Evidently 
many structural possibilities will have to be treated for a full description of ferric ion specia- 
tion in water, eventually considering higher aggregates and anharmonic vibrational motions 
of the strongly interacting water molecules and other ligands. These issues are probably 
best pursued through development of a molecular mechanics force field to screen structures 
rapidly, reserving electronic structure calculations for verification of the important struc- 
tures found and refinement of the force field. The energetic ordering found here for isomers 
of Fe(H 2 0) 4 (OH) 2 + is cis (lowest), outer sphere, and trans (highest), but all these energies 
are within about 2 kcal/mol. Molecular mechanics force fields might be reparameterized 
to account for these results 0,0,0 • Identification of prominent isomers should simplify and 
improve the modeling of the temperature variations of thermodynamic properties. It de- 
serves emphasis also that substantial contributions to these reaction free energies, and to 
the accumulated uncertainties, are associated with the water partners in these reactions, i.e. 
the free energies associated with solvation of H2O and H30 + . A similar comment would 
apply for other common oxy-acids and ligands in water, e.g. carbonate, nitrate, sulfate, 
and phosphate. Molecular descriptions of metal ion speciation will be incomplete without 
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accurate molecular characterizations of these species in aqueous solution. 
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APPENDIX 

Here we specify with greater care some of the statistical thermodynamic considerations 
relevant to solvation free energies obtained from cluster calculations of the present variety 
| |13| |. We note that these issues are of minor importance for the hydrolysis reactions that 
are the focus of this paper. However, these considerations become more important for the 
absolute free energy reported for the aqueous ferric ion and for the free energy of dissociation 
of the outer sphere complex. 

Utilizing calculations on clusters 

We first specify more fully and explain the procedure from computing the absolute free 
energy of the ferric ion given in Table IV on the basis of cluster results obtained from 
density functional theory and the dielectric continuum estimate of solvation contributions. 
The result of the development here is a formula for the chemical potential of the ferric ion: 

^ e3+ = RT H m{£-^))J - 6fiw > (7) 

Here pF e 3+ the number density of ferric ions in the solution, V is the volume of the sys- 
tem, q M is the partition function of the isolated hexaaquoferric complex (M=Fe(H 2 0) 6 3+ ), 
—RT\n{(e~ AU / RT ))o ! M is the solvation free energy of the complex, and is the chemical 
potential of the water. The conclusion drawn from this formula is that the desired absolute 
free energy of the hydrated ferric ion is obtained from the chemical potential change of the 
reaction Fe 3+ +6H 2 — > Fe(H 2 0)6 3+ except with the additional provisions that the chemical 
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potentials should be evaluated on a fully molecular basis at the water concentration of in- 
terest and no solvation contribution should be included for the atomic ion. The subsequent 
section notes how the isolated molecule partition functions should be used to obtain the 
chemical potentials at the required concentrations. 

As suggested above, the fundamental issue is the use of the cluster calculation Fe(H 2 0) 6 3+ 
to obtain the chemical potential, /iF e 3+, of the hydrated ferric ion, Fe 3+ . It is clear on physical 
grounds that such an approach should be advantageous when the identification of a relevant 
cluster is physically obvious and when the inner shell of the cluster requires a specialized 
treatment. What happens when the relevant cluster is not so obvious? What about cases 
when more than one cluster should be considered? How might this approach be justified 
more fully and how might the calculations be improved? 

The statistical mechanical topic underlying the considerations here is that of association 
equilibrium f45], 46l , f47 | and is often associated with considerations of 'physical clusters.' A 



suitable clustering definition [0,fEJ] is required for these discussions to be explicit and heavier 
statistical mechanical formalisms can be deployed P5]|iS|Ji7| . However, the treatment here 
aims for maximal simplicity. This argument is an adaptation of the potential distribution 
theorem p]|]. 

In order to involve information on clusters, we express the density of interest in terms 
of cluster concentrations. Thus, if the Fe 3+ ion appears only once in each cluster, i.e. if 
mononuclear clusters need be considered, then we would write 

PFe*+ = ^2pm (8) 

M 

where M identifies a molecular cluster considered and the sum is over all molecular clusters 



that can form. A satisfactory clustering definition {|5],fJ9| insures that each ferric ion can be 
assigned to only one molecular cluster, e.g. that a molecular cluster with one ferric ion and 
six water molecules is not counted as six clusters of a ferric ion with five water molecules. 
The calculations above assumed M = Fe(H20)6 3+ , and that was it. In the more general 
case that not all clusters are mononuclear, Eq. § would involve the obvious stoichiometric 

15 



coefficients. The concentrations pu are obtained from 

Pm = z Fe s + z^(q M /V) ((e- AU/RT )) QM • (9) 

i%m is the number of water molecules in the cluster of type M, (six in the example carried 
long here); z Fe 3+ and zw are the activity of the ferric ion and the water, respectively; i.e. 
z 7 = e^~<l RT ] qM = Qm(T) is a conventionally defined canonical partition function for the 
cluster of type M p5| , |46|^7| , |49[| . The indicated average utilizes the thermal distribution of 
cluster and solvent under the conditions that there is no interaction between them. AU is the 
potential energy of interaction between the cluster and the solvent. In the example carried 
along here we do not pay attention to any counter-ions since those issues are tangential to 
the current considerations. 

Eq. [| is most conveniently derived by considering a grand ensemble. Suppose we have a 
definite clustering criterion: a cluster of a ferric ion and tim water molecules is formed when 
exactly um water molecules are within a specified distance d of a ferric ion. In the example 
we have been carrying along, water molecules with oxygens within about 2.2 A of a ferric 
ion are in chemical interaction with the ferric ion. It would be natural to specify d < 2.2 A 
for clustered ferric- water (oxygen) distances. The average number < Nm > of such clusters 
is composed as 

E(z Fe 3+,z w ,T,V) < N M > = zp^+z^ 1 (10) 

£ N Fe3+ z^r- l ( Nw )z^- nM Q(N,V,T\n M + 1) 

N Fe3+ >l,N w >n M \nu J 

Here E(z Fe 3+, zw, T, V) is the grand canonical partition function; N Fe 3+ is the number of 
ferric ions in the systems and Nw is the number of water molecules; Q(N, V,T\um + 1) is 
the canonical ensemble partition function with one specific ferric ion and n^j specific water 
molecules constrained to be clustered. The binomial coefficient provides the number of 
riM-tuples of water molecules that can be selected from Nw water molecules. Because of the 
particle number factors in the summand the partition function there can also be considered 
to be the partition function for N — 1 ferric ions and Nw — um water molecules but with 
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an extra hm + 1 objects that constitute the cluster of interest. A reasonable distribution 
of those Hm + 1 extraneous objects is the distribution they would have in an ideal gas 
phase; the Boltzmann factor for that distribution appears already in the integrand of the 
Q(N, V, T\um + 1) and the normalizing denominator for that distribution is um^(1m(T). The 
acquired factor nj^- cancels the remaining part of the combinatorial Adjusting the 

dummy summation variables them leads to Eq. |9|. 

If we were to identify an activity for an M-cluster as % = Zp e 3 + z w I we would obtain 
from Eq. || the general statistical thermodynamic formulae of Reference |44[ . 



A virtue of the derivation of Eq. || sketched here is that the primordial activities z Fe 3+ and 
zy/ are clear from the beginning. This helps in the present circumstance where concentrations 
and chemical potentials of many other species and combinations will be of interest also. 

Combining our preceding results, we obtain 



E^WVl((e- A ^n\ . (11) 



PFe 3 + 

— = ^-^/"n\ e ' // ,M 

This is a reexpression according to clusters of a basic result known both within the context of 



the potential distribution theorem |5(J and diagrammatic (mathematical) cluster expansions. 



For the latter context see Eq. 2.7 of Reference JIT]]. Rearranging, we obtain the desired 
chemical potential 

n M 

^ + = -^lnE(^)((e-/-)) ,J. (12) 

This is the result that was sought. If higher-order clusters had been considered with a 
suitable clustering definition, the final result would have involved a more general polynomial 
of z Fe 3+; higher powers of z Fe 3+ would appear in the sums that would replace Eq. || because of 
the presence of higher powers of z Fe 3+ in some instances of Eq. |9|. A virtue of Eq. [12] result 
is that the thermodynamic activity of the water appears explicitly and that contribution 
may be included in a variety of convenient ways, perhaps utilizing experimental results. 
Furthermore, we see that there is no question whether a particular standard state for the 
water is relevant or, for example, whether only the excess part of the chemical potential of 
the water is required. 
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Notice that if the cluster definition had been restrictive enough that only the atomic 
M=Fe 3+ were present in appreciable concentration, i.e. 71^=0 for all clusters that need 
be considered, then Eq. [12] produces the previously known general answer for the chemical 
potential of a ferric ion in solution P^J50|] : 



fi Fe3+ /RT = \n[p Fe3+ V/q Fe 3 + ] - In (e- AU/RT ) Q (13) 

It would be natural to choose the electronic energy of the atomic ferric ion as the zero of 
energy and to anticipate that the degeneracy of the electronic degrees of freedom will be phys- 
ically irrelevant. Then it would be sufficient to put V/qp e 3+ = A|, e3 + = (h 2 /2%mF e 3+RT) 3 / 2 , 
the cube of the deBroglie wavelength of the ferric ion. mp e 3+ is the mass of the ion, and h is 
the Planck constant. In any case, we define the absolute free energy of the hydrated ferric 
ion as the second term on the right A/i^ e 3+ = —RT In (^e~ AU / RT ^. In view of Eq. [7] we then 
have 

W = RT H^~] ~ RT In ((e- AU/RT )) 0M ~ ^w, (14) 

where M=Fe(H 2 0)6 3+ . This Eq. [TJjis the formula that was used. 

A generalization of interest is the case that the solvent contains more than one species 
that may complex with a specific metal ion. For example, suppose that ammonia may be 
present in addition to water, that mixed complexes may form with ferric ion, and that these 
complexes have been studied as clusters in the same way that the hexaaquoferric ion complex 
was studied above. Then the result Eq. |l| is straightforwardly generalized by including the 
proper combinations of activities of the additional ligands possible. 

Standard State Modifications 

Many ab initio electronic structure packages, such as the Gaussian94 package used here, 
can produce molecular (or cluster) partition functions qM = Qm(T) and on this basis free 
energies for the species and reactions considered. It should be emphasized that these are 
typically applicable to a hypothetical ideal gas at concentrations corresponding to pressure 



p=l atm, see Table III. Thus, for example, these results determine the ideal chemical 
potential Hw = RTln[p w V/qw(T)]. However, because of the choice of hypothetical p=l atm 
ideal gas standard state, those results use pw = p/RT with p=l atm. In Gaussian94, for 
example, this logarithmic concentration dependence is considered a translational entropy 
contribution and, therefore, the entropic contribution to the reaction free energy of the first 
reaction in Table III is substantial and unfavorable. Because of the dilute reaction medium 
associated with p=l atm, the water molecules written as reactants in this reaction have more 
freedom before complexation than they do after. This is an entirely expected physical effect 
but inappropriate in solution. To obtain results applicable at the concentration of liquid 
water we determine that pressure parameter p = pwRT from the experimental density of 
liquid water pw= 997.02 kg/m 3 . The required value is p=1354 atm as indicated in Table IV. 
When this value is utilized in the expression for the translational entropy contributions, the 
translational entropy penalty of Table III for the first reaction is about half recovered. 
Because the second through fifth reactions of Table III were written to have the same 
molecule numbers for reactions and products this translational entropy contribution is very 
minor in those cases. 

We note also that the experimental tabulation of Reference |Ttl suggested as an exper- 



imental result for the absolute free energy of the hydrated ferric ion (-1037 kcal/mol) of 
Table IV requires an adjustment of about 1.9 kcal/mol |f24| , [42| for similar reasons. This 
adjustment is insignificant for that property with the methods used here. 
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FIGURES 



Figure 1. Computed structure of the Fe(H 2 0)6 complex. 

Figure 2. Computed structure of the Fe(H 2 0) 5 OH 2+ complex. 

Figure 3. Computed structures of the Fe(H 2 0) 4 OH 2+ complexes; cis (upper), outer 
sphere (middle), trans (lower), in order from lowest to highest energy. 

Figure 4. Computed structure of the penta-coordinate Fe(H20)3(OH) 2 f complex. This 
should be compared to the middle structure of Figure 3. 
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TABLES 



TABLE I: Electronic energy (E/[au]), Zero-point energy (ZPE/[kcal/mole]), and excess 
chemical potential (//*/[kcal/mol]) with the B3LYP and dielectric continuum solvation ap- 
proximations. Dielectric radii for all atoms were taken from Ref. JHJ except Rp e =2.08A. 





E 


ZPE 


* 

(J- 


H 


-0.497885 


- 


- 


Fe 3 + 


-1261.590210 


- 


- 


H 2 


-76.434010 


13.3 


-8.3 


H 3 0+ 


-76.707704 


21.5 


-96 


OH 


-75.739116 


5.2 




OH 


-75.802535 


4.5 


-108 


Fe(H 2 0) 6 3 + 


-1721.261664 


94.3 


-441 


Fe(H 2 0) 5 OH 2 + 


-1721.216989 


85.6 


-203 


cis Fe(H 2 0) 4 (OH) 2 + 


-1720.987174 


79.1 


-71 


trans Fe(H 2 0) 4 (OH) 2 + 


-1720.983894 


78.3 


-74 


outer Fe(H 2 0) 4 (OH) 2 + 


-1720.985680 


79.0 


-73 


Fe(H 2 0) 3 (OH)+ 


-1644.536430 


62.8 


-69 
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TABLE II: Gas-phase thermochemistries (Ai?o,kcal/mole) in the B3LYP approximation. 
Note that these include zero-point energy. The ionization potential (IP) of hydrogen and 
electron affinity (EA) of the hydroxyl radical are in eV. 

AE Expt. IIJ55] 



H 2 -> H + OH 115.5 118.0 

H 2 -> H+ + OH" 387.5 389.3 

H 2 + H+ -> H 3 0+ -163.5 -165.1 

2H 2 -> H 3 0+ + OH" 223.9 224.0 

IP(H) 13.55 13.6 

EA(OH) 1.76 1.83 



TABLE III: Ideal gas thermochemistries (kcal/mole) at T=298.15 K and p=l atm, uti- 
lizing the B3LYP approximation. The first column gives the electronic energy contribution, 
the second includes the zero-point energy, the third evaluates the energy at 298. 15K, the 



fourth gives the enthalpy at 298. 15K, and the fifth column is the Gibbs free energy. 





AE e 


AE 


A£ 298 






Fe 3 ++6H 2 -> Fe(H 2 0) 6 3+ 


-669.8 


-655.2 


-657.2 


-660.8 


-601.6 


Fe(H 2 0) 6 3+ + H 2 -» Fe(H 2 0) 5 OH 2 + + H3O+ 


-143.7 


-144.2 


-143.2 


-143.2 


-148.4 


Fe(H 2 0) 5 OH 2 + + H 2 -> cis Fe(H 2 0) 4 (OH) 2 + + H3O+ 


-27.5 


-25.8 


-26.6 


-26.6 


-25.7 


Fe(H 2 0) 5 OH 2 + +H 2 -> trans Fe(H 2 0) 4 (OH) 2 + + H3O+ 


-25.5 


-24.5 


-25.0 


-25.0 


-24.9 


Fe(H 2 0) 5 OH 2 ++H 2 -> outer Fe(H 2 0) 4 (OH) 2 + + H3O+ 


-26.6 


-25.0 


-25.5 


-25.5 


-25.9 


outer Fe(H 2 0) 5 OH 2 + -> Fe(H 2 0)3(OH) 2 + + H 2 


9.6 


6.7 


6.7 


7.3 


-1.8 
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TABLE IV: Aqueous reaction thermochemistries (kcal/mole) at T=298.15 K utilizing 
the B3LYP approximation. The first column gives the gas-phase free energy from Table III, 
the second reports the net free energy after adjustment for the alternative concentration in 
the liquid ("p=1354 atm" for the water species), and the third column gives the solvation 
increment to the free energy , i.e. that contribution due to solute-solvent interactions. 
Experimental results are given in the final column. The solution phase results of the first 
row are the "absolute free energies of the hydrated ferric ion." See the appendix, Eq. (0) 



for interpretation. 







^^298 


AG 2 98 




AG 298 + n* 


Expt. [0j55 


Fe 3 ++6H 2 -» Fe(H 2 0) 6 3+ 




-601.6 


-628.9 


-391.2 


-1020.1 


-1019,-1037 


Fe(H 2 0) 6 3 + + H 2 -» Fc(H 2 0) 5 OH 2 + + H3O+ 




-148.4 


-148.7 


150.3 


1.6 


3 


Fe(H 2 0) 5 OH 2 + + H 2 -» cis Fe(H 2 0) 4 (OH) 2 + + 


H3O+ 


-25.7 


-26.0 


44.3 


18.3 


5 


Fe(H 2 0) 5 OH 2 + +H 2 -» trans Fe(H 2 0) 4 (OH) 2 + 


+ H3O+ 


-24.9 


-25.2 


41.3 


16.1 




Fe(H 2 0) 5 OH 2 ++H 2 -> outer Fe(H 2 0) 4 (OH) 2 + + 


H3O+ 


-25.9 


-26.1 


42.3 


16.3 




outer Fe(H 2 0) 5 OH 2 + -» Fe(H 2 0) 3 (OH) 2 + + H 2 




-1.8 


+2.7 


-4.3 


-1.6 





TABLE V: Selected geometrical parameters for species in Fe(H20)6 hydrolysis reactions 
from B3LYP calculations. 





R(Fe-OH 2 )A 


R(Fe-OH)A 


#(OH-Fe-OH) 


0tfft(H 2 O)0 


Fe(H 2 0)3+ 


2.060 






0(6) 


Fe(H 2 0) 5 OH 2 + 


2.103 - 2.150 


1.760 




0(3),30(2) 


cis Fe(H 2 0) 4 (OH)^+ 


2.172 - 2.296 


1.820 - 1.847 


108 


36,40,44,51 


trans Fe(H 2 0) 4 (OH)^ + 


2.168 - 2.194 


1.851 


180 


38(2),44(2) 


outer Fe(H 2 0) 4 (OH)^+ 


2.108 - 2.196 


1.811 


114 


18,44,44 


Fe(H 2 0) 3 (OH) 2 


2.102 - 2.185 


1.803 - 1.807 


128 


22,40,40 



1 Defined as 0-Hi-H 2 -Fe dihedral angle. 
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TABLE VI: Atomic charges from ChelpG analysis of the electrostatic potential. 



Fe(H 2 0) 6 Fe(H 2 0) 5 OH Fe(H 2 0) 4 (OH) 2 Fe(H 2 0) 4 (OH) 2 Fe(H 2 0) 4 (OH) 2 Fe(H 2 0) 3 (OH) 2 









cis 


trans 


outer 




Fe 


+2.15 


+1.62 


+1.14 


+ 1.54 


+1.24 


+0.92 


02 


-1.02 


-0.88 


-0.69 


-0.85 


-0.69 


-0.63 


H8 


+0.58 


+0.50 


+0.43 


+0.47 


+0.41 


+0.41 


H9 


+0.58 


+0.52 


+0.42 


+0.48 


+0.42 


+0.42 


03 


-1.02 


-0.88 


-0.68 


-0.89 


-0.69 


-0.62 


H10 


+0.58 


+0.50 


+0.43 


+0.49 


+0.41 


+0.41 


Hll 


+0.58 


+0.52 


+0.42 


+0.48 


+0.42 


+0.41 


04 


-1.02 


-0.78 


-0.62 


-0.86 


-0.66 


-0.48 


H12 


+0.58 


+0.47 


+0.40 


+0.45 


+0.42 


+0.37 


H13 


+0.58 


+0.47 


+0.39 


+0.48 


+0.42 


+0.37 


05 


-1.02 


-0.77 


-0.71 


-0.83 


-0.88 


— 


H14 


+0.58 


+0.47 


+0.43 


+0.45 


+0.45 


— 


H15 


+0.58 


+0.47 


+0.42 


+0.48 


+0.45 


— 


06 


-1.02 


-0.80 


-0.82 


-0.90 


-0.77 


-0.67 


H16 


+0.58 


+0.47 


+0.44 


+0.45 


+0.40 


+0.39 


H17 


+0.58 


+0.41 










07 


-1.02 


-0.97 


-0.82 


-0.89 


-0.75 


-0.68 


H18 


+0.58 


+0.60 


+0.44 


+0.44 


+0.38 


+0.39 


H19 


+0.58 













28 




Figure 1 




Figure 2 




Figure 3 




Figure 4 



